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Employing a medium-dependent quantum optics formalism and a Green function solution of 
Maxwell's equations, we study the enhanced spontaneous emission factors (Purcell factors) and 
Lamb shifts from a quantum dot or atom near the surface of a slow-light metamaterial waveguide. 
Purcell factors of approximately 250 and 100 are found at optical frequencies for p— polarized and 
s— polarized dipoles respectively placed 28 nm (0.02 Ao) above the slab surface, including a realistic 
metamaterial loss factor of r y/2ir = 2THz. For smaller loss values, we demonstrate that the slow- 
light regime of odd metamaterial waveguide propagation modes can be observed and related to 
distinct resonances in the Purcell factors. Correspondingly, we predict unusually large and rich 
Lamb shifts of approximately —1 GHz to —6 GHz for a dipole moment of 50 Debye. We also make 
a direct calculation of the far field emission spectrum which contains direct measurable access to 
these enhanced Purcell factors and Lamb shifts. 



I. INTRODUCTION 



Early in 1968, Veselago predicted that a planar slab 
of negative-index material (NIM), which possesses both 
negative permittivity s and negative permeability ^, 
could refocus electromagnetic waves While an in- 
teresting prediction, these so-called metamatcrials did 
not receive much attention in optics research until Sir 
John Pendry showed that NIMs could be used as a "su- 
perlens", which could overcome the diffraction limit of 
conventional imaging system Subsequently, a host 
of applications has been proposed, ranging from designs 
for optical cloaks to hide objects 0] , through to schemes 
that completely stop light Q. Although most of these 
schemes are idealized, and suffer in the presence of meta- 
material loss, they have nevertheless motivated signifi- 
cant experimental progress. For example, Schurig et al. 
experimentally demonstrated some cloaking features by 
split-ring resonators [H ; and recent achievements in fabri- 
cation have facilitated the realization of negative indices 
at communication wavelength [||, with some extension 
to quasi-3D structures also reported 0]. 

While it is certainly becoming established that meta- 
materials posses some remarkable classical optical prop- 
erties, less well studied arc their quantum optical proper- 
ties, such as what happens to the spontaneous emission 
of an embedded atom or quantum dot. In 1946, Purcell 
pointed out that due to the spatial variation of the local 
photon density of states (LDOS), the spontaneous emis- 
sion rate in a cavity can be enhanced or suppressed de- 
pending upon the distance between the mirrors Q . The 
modification of spontaneous emission due to inhomogc- 
neous structures is a large field in its own right, leading 
to applications in quantum optical technology [9j . In the 
domain of metamaterials, Kastel and coworkers inves- 
tigated the spontaneous emission of an atom placed in 
front of a mirror with a layer of NIM [l(| ; this study was 
motivated by the perfect lens prediction of a vanishing 
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optical path length between the focal points, leading to 
the peculiar property that the evanescent waves emerg- 
ing from the source are exactly reproduced; consequently 
they found that the spontaneous emission can be com- 
pletely suppressed. This prediction does not account for 
the essential inclusion of metamaterial loss or absorption; 
yet, it is well known that metamaterials must be disper- 
sive and absorptive to satisfy the fundamental principle 
of causality and the Kramers-Kronig relation [111 ]. Not 
surprisingly, when absorption is necessarily taken into 
account, then the predicted properties of an ideal loss- 
less metamaterial can qualitatively change. For exam- 
ple, the superlens and the invisible cloak are never per- 
fect [HI, [l3[, and slow light modes can never be really 
stopped and are usually impractically lossy Simi- 
larly, it is expected that absorption will have an impor- 
tant influence on the modification of spontaneous emis- 
sion 0]. In this regard, Xu et al. have extended the 
works of Kastel et al. to one dimensional right-handed 
and left-handed material layers and find that nonradia- 
tive decay and instantaneous radiative decay will cer- 
tainly weaken the predicted inhibition of spontaneous 
emission (iH [l7|. 

Some of the first theories to treat quantum electro- 
dynamics near a interface were introduced around 1984 
by Wylie and Sipe Ed, where, using Green function 
techniques, they showed that the scattered field can be 
expressed in terms of the appropriate Fresnel suscepti- 
bilities. Using such methods, it is now well known that 
the photonic LDOS can be increased near a metallic sur- 
face, e.g. see Ref. [2(|, whereby the p— polarized dipole 
couples to a transverse magnetic (TM) surface plasmon 
polariton (SPP). Typically such resonances are far from 
the optical frequency domain and they are restricted to 
TM polarization; in addition, the emission is dominated 
by quenching or non-radiative decay. Even in the pres- 
ence of gratings, enhanced emission via SPPs is not very 
practical [2lj . However, the rich waveguide properties 
of metamaterials have quite different polarization depen- 
dences and mode structures than SSPs at a metal surface; 
for example, they can support slow-light, bound propa- 
gation modes and transverse electric (TE or s— polarized) 
SPPs. It is therefore of fundamental interest to explore 



2 



the quantum optical aspects of these novel waveguides. 

Enhanced emission at the surface of both metals and 
metamaterial slabs was studied in 2004 by Ruppin and 
Martin (22j. They noted that resonance peaks due to 
s— polarized surface modes and waveguide modes can ap- 
pear for the metamaterial case, although they did not 
discuss the origin of the waveguide peaks. Similar find- 
ings were found by Xu et al. HH, but the role of loss 
was not explored in detail but rather it was treated as 
a perturbation and assumed to lead to only dissipation; 
such an assumption is highly suspect in a NIM waveg- 
uide, since the entire modal characteristics of the struc- 
ture depend intimately upon the material loss and dis- 
persion profiles [l4| . Very recently, spontaneous emission 
enhancements in NIM materials and interesting quantum 
interference effects have been reported by Li et al. (27| . 
although unrealistically small losses were typically as- 
sumed, and again the physics behind the enhancement 
factors was not made clear. In all of these works, there 
has been no quantitative connection made to the complex 
band structure or to the far-field (and thus measurable) 
spontaneous emission spectra or dipolar frequency shifts 
(Lamb shifts). 

In this work, the modified spontaneous emission of a 
quantum dot or atom (single photon emitter) situated 
above a slow-light metamaterial waveguide is investi- 
gated in detail by employing a medium-dependent Green 
function theory and comparing with the lossy guided 
waveguide modes. We compute the PF as well as the 
spontaneous emission spectrum in the far field by devel- 
oping and using a rigorous quantum optics theory We 
stress that the recent prediction of completely stopped 
waveguide modes in a metamaterial waveguide [H would 
lead to an infinite PF, but as reported by Reza et al. [bH ]. 
the properties of the slow-light modes are significantly 
changed in the presence of loss; thus we also investigate 
the dependence on loss in some detail. We show that 
the emission properties of a photon emitter can act as a 
probe for below light-line propagation mode characteris- 
tics, showing measurable enhanced radiative broadening 
and quantum Lamb shifts. The Lamb shifts are found to 
be extremely rich and pronounced. We also compare and 
contrast these NIM quantum optical features with well 
known results for metallic surfaces. 

Our paper is organized as follows. In Sect, [ill we intro- 
duce the NIM waveguide structure of interest and com- 
pute and discuss the band structure for both s (TE) and 
p (TM) polarization. In Sect. IIIH we present a rigor- 
ous theory for calculating the Purccll factor (PF) , Lamb 
shift, and emitted spectrum from a single photon emit- 
ter. From this theory, we derive an explicit and analytical 
solution to the emitted field at any spatial location us- 
ing a full non-Markovian theory which is valid for any 
general media (lossy and inhomogeneous). In Sect. IIV1 
we discuss a general technique for computing the multi- 
layered Green function using a stratified medium tech- 
nique of Paulus et al. [28| , and formally separate the total 
Green function into a homogeneous and scattered part to 
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FIG. 1: Schematic diagram of the system being investigated. 
The green dot at r = refers to the atom or quantum dot 
that is decaying radiatively (with a rate T), at a distance zo 
above the slab of thickness h. The vertical position of the 
dot, z d = z . 



properly obtain the photonic Lamb shift. In Sect.fVj we 
present calculations for the Purccll effect and Lamb shift , 
as well as the spontaneous emission spectra emitted into 
the far field. Finally, we give our conclusions in Sect. I VII 



II. METAMATERIAL SLAB WAVEGUIDES: 
COMPLEX BAND STRUCTURE AND 
SLOW-LIGHT PROPAGATION MODES 

In the following sections, we wish to calculate the spon- 
taneous emission from dipoles d = de a , where a = z or 
a = x(y). Because the spontaneous emission is strongly 
affected by the metamaterial waveguide modes, in this 
section we present the results of calculations of the dis- 
persion of both TE and TM propagation modes. The 
schematic diagram of the system under study is shown 
in Fig. [TJ The negative index slab is surrounded by air 
and assumed infinite (or much larger than a wavelength) 
in the x and y directions. The thickness of the slab is 
h = 280 nm. In view of the importance of dispersion and 
absorption, we take both into account from the begin- 
ning, which ensures causality and thus avoids unphysical 
results. The dispersion is introduced via the Lorentz and 
Drudc models 1291] : 



H{uj) = 1 + 



pm 



e(w) = 1 



LU 2 + ijLU 



(1) 

(2) 



where uj pm and ui pe are the magnetic and electric plas- 
mon frequencies, and luq is the atomic resonance fre- 
quency. In what follow, we are interested in waveguides 
with guided modes at optical frequencies and thus take 
iv /2n = 189.4 THz, iv pm /2n = 165.4 THz, uj pe /2ir = 
490 THz. 

To solve the complex band structure, one can work 
with a complex wave vector (/?) and a real frequency (u) , 
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FIG. 2: Dispersion curves of the lossy metamaterial waveg- 
uide for the first few lower-order modes using the complex- 
es approach for -y/2-K = 2 THz (7/wo ~ 0.01). The red 
and blue curves show the surface plasmon polariton modes 
and bound propagation modes, respectively. The solid and 
dashed curves represent the even and odd modes, respecitvely. 
The solid thick black curves display the vacuum light-line and 
metamaterial light-line, while the horizontal dashed lines in- 
dicate the atomic (wo) and TE plasmon resonances (w^). (a) 
Re(aj) versus j3 for TM modes. The modes become more dense 
near the resonance frequency Wo, and form a continuum, (b) 
Im(w) versus j3 for TM modes. (c)(d) The same as (a)(b) but 
for TE modes. 



or, alternatively, a complex— to and a real— (3. The for- 
mer is perhaps more appropriate for modeling plane wave 
excitation, while the latter is better suited for a broad- 
band excitation response that is invariant in z. Neither 
of these approaches constitutes a complete connection to 
a broadband dipole response, and thus we will show both 
solutions, and briefly discuss their main features. The de- 
tails of these two approaches for modelling metamaterial 
waveguide properties will be presented elsewhere. 

The complex dispersion curves of the aforementioned 
metamaterial waveguide for both TE and TM modes are 
shown in Fig.f2]and Fig. O using complex-w and complex- 
f3 approaches, respectively. In the complex—/? approach, 
we only show the Rc[/3] solution (Fig. [3]), as the Im[/3] 
simpl y d emonstrates the large losses in the regime of slow 
light [14| . These curves come from the complex solution 
to the transcendental dispersion equation derived from 
the Maxwell's equations and the guidance conditions [3fj| . 

In Fig. we show the first few TM (TE) propaga- 
tion modes: TM2-TM6 (TE2-TM6), as well as TE0-TE1 
(TM0-TE1) SPP modes (red curves); a wide range of 
frequencies is displayed from from 180 to 240 THz. We 
stress that the TE-polarized SPPs are unique to meta- 
materials and the properties of these SSP modes can be 
engineered to have a resonance in the optical frequency 
regime, near the propagation modes. For the TM case, 
there is also a higher-lying SPP resonance, which we do 
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FIG. 3: Dispersion curves of the lossy metamaterial waveg- 
uide for the first few lower-order modes in using the complex-/? 
approach, where only Re[/3] is shown, (a) TM modes, (b) TE 
modes. The thin solid black curves represent vacuum light- 
line and metamaterial light-line. Note that formally labelling 
the curves is more difficult than in Fig. 2 due to the merging 
of the different curves (e.g. the TM1 and TM5 curves near 
220 THz). 



not show as it is far outside the frequency of interest 
(w*™ - 346 THz or specifically u pe /V2); the TM0 and 
TM1 modes show the start of the SPP modes just below 
the air light line. 

The properties of these metamaterial waveguide modes 
shown in Fig. [2] are considerably different from those in a 
conventional dielectric waveguide [3l| . The most impor- 
tant difference is that backward and forward propagat- 
ing modes exist and, in a lossless metamaterial waveg- 
uide, even stopped-light modes (where the slope goes 
to zero) are supported. However, the necessary inclu- 
sion of intrinsic loss in a metamaterial slab dramatically 
changes the dispersion curves near slow-light regions (see 
also Ref. [U[). Thus one must include material losses to 
have any confidence in the results and predictions. In 
the complex— u> solutions (Fig. 2), although the slope and 
thus the group velocity v g is zero at some points, the use 
of the group velocity as a measure of energy transport 
is not meaningful [14| . and there is a finite imaginary 
part of frequency. We note, however, that at the points 
where the slope is zero in Fig. HJa) and (c), the density 
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of modes is large and as we shall see, this can lead to a 
large Purcell factor at the relevant frequencies. We note 
in particular that because the slope of all of the different 
mode dispersions tends to zero as /? goes to infinity, there 
is a large enhancement in the density of states at lo = lo$. 

In Fig. [3] we show the dispersion curves using the 
complex—/? approach. In this approach, the impact of 
material losses on the dispersion curves are much more 
drastic and these curves would appear to have little to do 
with the curves in Fig. [21 However, despite their differ- 
ences, there is a strong correspondence between the two 
sets of curves. In a (fictitious) lossless waveguide, the 
two sets of curves would be identical apart from some 
leaky modes that carry no energy. In the complex—/? 
approach, two degenerate leaky modes start at the zero- 
slope point of the non-SSP propagation modes and move 
to higher frequencies. In Fig. [Jl these leaky modes split 
and merge seamlessly into the propagation mode dis- 
persion, creating the split-curve structure that is seen, 
e.g. for the TE4 mode near lo/2tt = 201 THz. Although 
the slope no longer goes to zero in the complex—/? ap- 
proach at the point where the propagation mode splits 
into the two leaky modes, the energy velocity (which is 
the correct measure of energy transport in a lossy sys- 
tem) is quite small but nonzero. Another important dif- 
ference between the modes in the complex—/? approach 
and the complex-w approach is that for the modes in 
the complex— a; approach, /? — > oo as lo — > lo$, while in 
the complex—/? approach, the modes bend back on them- 
selves at the atomic resonance frequency loq. Note that 
the group velocity, which is given by the slope of disper- 
sion curve, is infinite at this point. There is no paradox 
here because, as discussed above, the group velocity is 
ill-defined in a lossy waveguide structure and the en ergy 
velocity is a correct measure of the transport speed [321 ] . 
We have calculated the energy velocity and find that the 
minimum energy velocity occurs exactly at the resonance 
frequency u>o, but that it is never zero; for example, the 
energy velocity minimum for both TE3 and TM3 modes 
is found to be around 10~ 3 c to 10~ 4 c, where c is the 
vacuum speed of light. A similar effect occurs for TE 
modes near the plasmon resonance, where in the in the 
complex—/? picture, the SSP modes not only penetrate 
below the plasmon frequency, but they transform seam- 
lessly into the higher-order leaky modes. For example, 
the TE4 and TE5 modes merge into SSP modes. Inter- 
estingly, in Fig. 3(b), we also observe new resonances in 
the left branch of the leaky TE propagation modes, which 
couple to the TE SPP modes just below the bare SPP 
resonance. Thus, two resonances appear in this SPP fre- 
quency regime; evidently, this is not expected from the 
complex— u; perspective. 

Later, we will show that the resonant frequency and 
the slow light regimes of the propagation modes are ex- 
actly coincident to that of the LDOS peaks in the spon- 
taneous emission spectrum, and thus both complex— lo 
and complex—/? band structures are useful to gain in- 
sight into the origin of the LDOS peaks. The spatial 



symmetries of the even and odd modes will also prove 
to be very important; since the odd modes have a much 
larger field amplitude at the surface, they couple much 
more strongly to a quantum dot or atom near the NIM 
surface. 



III. QUANTUM THEORY OF SPONTANEOUS 
EMISSION IN METAMATERIAL 



A. Purcell factor (PF) and Lamb shift 

The PF is a measure of the spontaneous emission rate 
enhancement; it is defined as PF = T/To, where T (the 
Einstein A coefficient) is the spontaneous emission rate 
associated with population decay rate from an excited 
state to the ground state, and r is the spontaneous emis- 
sion rate in vacuum or a lossless homogeneous medium. 
We consider a system consisting of a quantum dot embed- 
ded in or near a general dispersive, absorptive, and inho- 
mogeneous medium. Employing a quantization scheme 
that rigorously satisfies the Kramers-Kronig relations, 
and using the electric dipolc approximation, an appro- 
priate Hamiltonian - following the works of Welsch and 
coworkers - can be written as [33, yj, ^al : 



tot) 



H =huj d a + a + ^2 J ' dr J dw;^f](r, coi) ■ f A (r, 



[a+d + a-d] -F(r d ), 



(3) 



where f|(r, lo{) and f\(r, LOi) are the continuum bosonic- 
field operators of the electric (A = e) and magnetic field 
(A = m) with eigenfrequency loi, a^' are the Pauli op- 
erators of the electron-hole pair (exciton), and iOd and 
d = ridd (d is assumed to be real) are the transition fre- 
quency and dipole moment of the dot, respectively. The 
field operator F is essentially the electric field operator 
augmented by the quantum dot polarization and can be 
expressed as F = T)/(eqe) + P/ (sos), where D is the dis- 
placement field (E = D/(e £) is the electric field away 
from the dipole), e is the complex relative permittivity 
and P is the polarization arising from the quantum dot 
dipole; this latter contribution is needed as it is the dis- 
placement field that should coup le to the dipole in the 
interaction Hamiltonian [l8|, [36[ |37|. Using the above 
formalism, we derive 



F (M) = hi— J dtoildr'G(r,r';ioi)- y/e^r',^) 



f e (r',cJi;£) H V x s/-m{r', loi) f m (r', u>r,t) 

e e(r) 

where the last term represents the polarization field from 
the dipole, Ej(r,LOi) and Kj(r, loi) are the imaginary parts 
of e(r,Wi) and l//i(r, loi) respectively, and e(r, loi) and 



fx(r,u!i) are the relative complex permittivity and perme- 
ability. The dyadic G(r, r';w;) is the electric-field Green 
function (GF) that describes the field response at r to 
an oscillating polarization dipole at r' as a function of 
frequency; the GF is defined from 



1 ujf 

V x — -V x — !re(u!i,r) 

H{u>i,r) c 2 



G(ry ;W |) = 



(5) 



where 1 is the unit tensor. 

From the above Hamiltonian, we derive the Hcisen- 
berg equations of motion for the time-dependent opera- 
tor equations as (t is implicit): 



da 

~~dT ~ 

df e (r,u>i) 
dt 



iujia +ih 1 d-F(r d ), 
= -iuJii e {r,LOi) 



(6) 



d-G*{r d ,r;uJi)[a- +a + ], (7) 



df m (r,wi) 
dt 



-iujif m (r,uji) 



d-[G'(r d ,T;u)i)xV][&-+& + ], 

(8) 



where [G*(r d ;r;w ; ) x V]„ = e jk id k G^(r d) T;ui), and we 
have used the one photon or weak excitation approxi- 
mation through a z F = — F. We can make a Laplace 
transform on the above set, defined through Oi(u>) = 
f °° e iuJt O l (t)dt, and obtain 



a (uj) 
a + (uj) 



i&-(t = 0) h~ x d ■ F(v dl ui) 

U) - U) d UJ - Ul d 

ia+(t = 0) rtifrj.w) 

UJ + UJ d 

f e °(r,w/;w) 



uj + UJ d 



(9) 
(10) 



' e/(r ' Wl) -d.G'(r <I> r;a„)- i[fr " (t,,) + fi+(w)] 



f m (T,u>i;uj) 



fm(r,w;;^) 
-ki(t,uji) c 



UJ -UJl 



(11) 



7r?ieo uji 
i[&-{uj) + a+(uj)} 



d [G*(r d ,r;w ; ) x V] 



UJ - UJl 



(12) 



where f represents a possible free field or homogeneous 
driving field in the absence of any quantum dot or atom. 

We next assume that the initial field is the vacuum 
field (i.e., f£(r,w z ;w) = 0), substitute Eqs. l|lU12[> into 



Eq. (4) , and make use of the relation (see [33| ) : 



ds{-Ki[G(r,s;uji) x V S ][V S x G*(s, r'; uji)}— + 

u i 

e/(s, wi)G(r, s, W| ) ■ G*(s, r'; wj)} = ImG(r, r'; uj x ), (13) 

Subsequently, we obtain an explicit solution for the dipole 
operators (and thus the polarization operator): 

o~ (uj) + & + (ui) = 

ia~{t = 0)(w + ui d ) + ia + (t = 0)(uj - ui d ) 



uj-uji- 2uj d d ■ [G(r d , r d] uj) + 5 -^-\ ■ d/% 



£0 



(14) 



which we can rewrite as 

<7~ (uj) + V + (uj) = 

ia~(t = 0)(uj + uj d ) + ia + (t = 0)(w - uj d ) 
ui d - uj 2 - 2uj d d ■ K(r d , r d ;uj) ■ d/he 



(15) 



where the new GF, K(r, r';uS) = G(r, t';u>) + S(r — 
r')/e(r). This is exactly the same form as the GF used in 
the formalism by Wubs et al. [38j , where the K function 
appears naturally when working with mode expansion 
techniques for lossless inhomogeneous media. The origin 
of the discrepancy between theories that use G or K, is 
because the correct interaction Hamiltonian should really 
contain a displacement field [TH, [36|, H3] , as has been ac- 
counted for in our Eq. (4). This subtlety becomes impor- 
tant, e.g., when deriving a Lippman-Schwingcr equation 
for the electric-fiel d op erator, which can only be achieved 
through use of K [331 • 

It is worth noting that the above equations are ob- 
tained with no Markov approximation, so they can be 
applied to both weak and strong coupling regimes of 
cavity-QED. In addition, we have made no rotating- wave 
approximations. In the weak to intermediate coupling 
regime, the decay rate of spontaneous emission V can be 
conveniently expressed via the photon GF through, 



r(r d) 



2d-Im[K(r d ,r d ;u;)] -d 

he 



(16) 



where for free space, Im[K vac (w)] = Im[G vac (cj)] = 
uj 3 /6ttc 3 , and so T = 2d 2 ui 3 /(he Q 6Trc 3 ). 

Within the dipole approximation, the above formal- 
ism is exact, and for lossless media, Eq. (|16[) can be 
reliably applied as soon as G is known, and one can 
exploit Im[G(r, r'; u>)] = Im[K(r, r';u>)], since e is real. 
In a previous paper dealing with lossless photonic crys- 
tals [39[, two of us adopted precisely K(r, r';uj), since 
it can be constructed in terms of the transverse modes. 
However, for lossy structures such as metals and meta- 
materials, Im[G(r, r; uj)] diverges [iF|, so we are forced 
to confront the immediate unphysical consequences of 
a dipole approximation. For both real and complex e, 
Rc[G(r, r; uj)] also diverges, which is well known. These 
GF divergences, as r — > r', are of course not physical 
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and arc simply a consequence of using the dipole approx- 
imation. Any finite size emitter, no matter how small, 
will have a finite PF and a finite Lamb shift [4l[. The 
usual procedure for a lossy homogenous structure is to 
either regularize the GF by introducing a high momen- 
tum cut off (4^ |. or to introduce a real or virtual cavity 
around a finite size emitter and analytically integrate the 
homogenous GF [43| . 

In the remainder of this paper, we will only concern 
ourselves with dipole emitters located in free space above 
a NIM waveguide; we will, however, revisit the problem 
of dipoles inside a NIM in future work, where one must 
carefully account for local field effects. For any inhomo- 
geneous structure such as a layered waveguide, a conve- 
nient approach to using the GF is to formally separate 
it into two parts, namely a homogeneous contribution 
Qhom ( wnose solution can be obtained analytically) and 
a scattering contribution G scatt . This approach, which is 
especially well-suited to dipoles in free space, is the ap- 
proach that we will adopt below. Using this separation, 
one can identify the non-divcrgent PF and Lamb shift 
solely from the scattered part of the GF. We obtain 

rsca t Vr 2d.Im[G scatt (r rf ,r d ^)]-d 

heo 

so that the Purccll factor, for a dot in free space, is 



PF{r d ) = 1 



fc (r d ^) 



r» 

Similarly, the Lamb shift is given by 



(18) 



a / x d-Re[G scatt (r d ,r )]-d 

Aoj(r d ,uj) = , (19) 

he 

where we have neglected the vacuum Lamb shift from 
the homogenous GF, since it can be thought to exist al- 
tready in the definition of u>d, and in any case it will be 
much smaller than any resonant frequency shifts com- 
ing from G scatt . While, in principle, one can apply mass 
renormalization techniques to also obtain the vacuum (or 
electronic) Lamb shift [44[, any observable shift will be 
related to - and completely dominated by - the photonic 
Lamb shift, and thus from G scatt . 



B. Spontaneous emission spectrum 

Next, we will obtain an exact expression for the emit- 
ted spectrum. From Eq. (4) and Eqs. (11-15), we obtain 
the analytical expression for the electric field operator, 
E(R,w), forR^r d : 

E(Rc) = / dun — ImG(R, r.^-d ^ 1 ^^' M 

= -G(R,r d ,u;)-d[a-(uj) + a + (io)}, (20) 

where -, — r = ttS(uii —uj)+iP( — - — ) (uii is the inte- 

gration variable) has been used. Note that the principal 



value of the integral cannot be neglected, otherwise only 
the imaginary part of GF in Eq. (20) is retained. This 
can be contrasted to the expression derived by Ochia 
et al. (45[, where the real part of the GF was omitted 
because they neglected the principal value of integral; 
however, this is unphysical and yields incorrect spectral 
shapes in general. 

The power spectrum of spontaneous emis- 
sion can be obtained from S(R.,uj) = 
Jo 00 /o° *i d< 2e^( t2 - tl )(E-(t 1 )E+(t 2 )), leading to 

S(t,u) = ((E(w))+E(w)). Using Eq. (20) and Eq. (15), 
one has, again for R^rj, 



d • G(R, r d ; + uj d )/e 



■ or 



2w d d- G scatt (r d ,r d ;w) ■ d/fie 



(21) 

This is in an identical form to the emission spectrum de- 
rived for a lossless material [3^] , showing that the electric- 
field spectrum at r depends on the two-space point GF, 
G(R, r d ]Lo), which describes radiative propagation from 
the dot position to the detector. All that remains to be 
done is obtain the GF, which we discuss next. 



IV. MULTI-LAYERED GREEN FUNCTION: 
PLANE WAVE EXPANSION TECHNIQUE 

The classical GF, G (r, r'; w), describes the response of 
a system at the position r to a polarization dipole located 
at r', so that the total electric field E (r, u) = G (r, r'; o)V 
vjx 1 ,ui) . In the case of a multilayer planar system [28|, 
|47| . when calculating the electric field in the same layer 
as the dipole, and as mentioned earlier, it is possible to 
formally write the GF in terms of a homogenous part and 
a scattered part. Formally one has 

G (r, r'; u) = G hom (r, r'; u) + G scatt (r, r'; u) . (22) 

Because we consider z and x(y) oriented dipoles sepa- 
rately, we only require the diagonal elements of the GF 
above the slab and the GF tensor elements can be greatly 
simplified. We take the source and field points to be 
r = (p,z) and r' = (p,z'), i.e. the transverse position, 
p = (x, y) of the dipole and observation points are equal. 
We will use the following notation to label the three- 
layered structure: region 1 is air, region 2 is metamatc- 
rial and region 3 is air. For the total GF, when both z 
and z' are in region 1, we have, for z < z' , 



^xx/yy (P, Z,Z ,W) = ^3 J 



— (r{ _ e lk ^( z+z ') +e - lk ^( z - z ')) 

-til f r P e ikiz(z+z') _|_ e -ifci« (*-«') 



(23) 
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G zz (p, z, z',uj) 



dk„ 



and with z > 



— — — r-d (z — z j H — 5 

e(r',w) 47tc 2 

ifcli(z+^') r g— ifcl*(z— z') 



A" 3 / 



(24) 



&a;x/ V y{P, z , z ,U) = = - / dk n k n 



j— (rf e ifc iz( z + 2 ') -|_ e *i 2 ( z_2 ')j 

^li gifci*(z+*') _ g^(z-z') 



G zz (p, z,z',w) 



fc3 



i\i (r', w) w 2 

W Jo ^ P hJ4 



(25) 



(26) 



Here, for s (TE) and p (TM) polarization, 



>/p) 



^{ s /p)^{ s /p)^{ s /p) a 2ik 2z h 



S s /p) I r 12 r 21 '23 



12 



1 _ r ( s /P) r ,( s /P)„2ife 2 ,/i 
1 r 21 r 23 e 



(27) 



and the wave vector, fc; = nuj/c in medium /. Cal- 
culating the z component of the wave vector as fc; z = 

1/2 

± (k 2 — k 2 ) when Re(fc;) > Re(k p ) where the positive 
sign is for positive index materials and the negative sign 
is for negative index materials. For Re(fc;) < Re(fc p ) we 

1 /2 

have h z = i {k 2 — kf) for both positive and negative 
index materials [33], |46| . The reflection and transmission 
coefficients arc 



t 



£jkiz ^ikjz 
£jkiz ~t~ &ikjz 

£jkiz H~ &ikj z 



(28) 



(29) 



From the above expressions it is seen that the scattered 
GF for any z in region 1 can be written as 



GIT* (p, *,*!» = 



ijj, (t',oj)uj 2 



Anc 2 



dkp kp 



r i ,\ 

hz c ik iz( z + z ) 

2k lz 



kur\ t 

~^2kJ 



%k\ z (z+z'} 



Gscatt / „ / \ 

zz \P,z,z,uj) 



i\x (r', bj) uj 1 
Aire 2 



dk n 



k\ z k^ 



P_ r P £ ik lz (z+z') 



(30) 



(31) 



2' 1 



We highlight again that the difference between Eqs. (1231 
and Eqs. (|30j|3ip is simply the homogeneous GF [471 ] . 



From a numerical perspective, the task is to solve Som- 
merfeld integrals. Such equations can be integrated in 
the lower half of the complex plane using the method 
described by Paulus et al. [28| for positive index ma- 
terials (where the poles are in the first and the third 
quadrant of the complex plane), or in the upper half of 
the complex plane for negative index materials (where 
the poles are in the second and the fourth quadrant of 
the complex plane); this method avoids any poles which 
may be near the real k p axis and improves numerical 
convergence, though it is unnecessary for large mate- 
rial loss. For our particular calculations, we integrated 
Eqs. (|23ti31[) using an adaptive Gauss-Kronrod quadra- 
ture which was verified to be well converged for a rela- 
tive tolerance of 10~ 4 . Specifically, the above equations 
were integrated along an elliptical path around the re- 
gion containing the bound and radiation mode contri- 
butions [28|, with the semi-major axis was 3 |Re (fe)) /2 
and the semi-minor axis was |Re (A?2 ) | /1000. After in- 
tegrating along the elliptical path, the equations were 
integrated into the evanescent region along the real k p 
axis. An additional advantage of this technique is that 
the integrand contributions from the bound and evanes- 
cent modes can be conveniently compared with the band 
structure, by examining the individual s— and p— polar- 
ized contributions as a function of k p for a given to, where 
it becomes obvious that the full GF solution requires both 
complex— uj and complex— (3 pictures. Since the GF ap- 
proach may be termed the complete answer, it is clear 
that the band structure approaches, either complex— a; 
or complex—/?, merely yield a limited sub-set solution 
about dipole coupling in these structures; having said 
that, both approaches (band structure and GF) offer a 
clear connection to the underlying physics of Purccll fac- 
tors and Lamb shifts. 



V. SPONTANEOUS EMISSION 
CALCULATIONS FOR A SLOW-LIGHT 
METAMATERIAL WAVEGUIDE 

A. Enhancement of the spontaneous emission rate 
(Purcell effect) 

The motivation behind investigating slow light waveg- 
uides in the context of enhanced spontaneous emission 
is that, quite generally, the relevant contribution to the 
LDOS from a lossless waveguide mode is inversely pro- 
portional to the group velocity [48], and so slow light 
modes may lead to significant PFs. In the field of planar 
photonic crystal waveguides, GF calculations [4!| and re- 
cent measurements [50j have obtained PFs greater than 
30 for group velocities that are about 40 times slower 
than c. This enhancement leads to an increase in the 
degree of light-matter interaction, and is important for 
fundamental processes such as nonlinear optics, and for 
applications such as single photon sources. The major 
difference with lossless photonic crystal waveguides and 
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NIM waveguides is that the N1M losses will likely mean 
that they are unlikely to find practical application as ef- 
ficient photon sources, since the photon emission is prob- 
ably dominated by non-radiative decay. Nevertheless, it 
is fundamentally interesting to calculate the emission en- 
hancements rates, and to connect these to a measurement 
that would allow direct access to this enhanced light- 
matter interaction regime. 

For our PF calculations, we first investigate the behav- 
ior of the spontaneous emission as a function of frequency. 
The GF is obtained directly using the multilayer GF tech- 
nique described above. Figure H] shows the spectral dis- 
tribution of the PF, when the emitter is placed above the 
slab (zd = z (see Fig. 1) = h/10 = 28 nm) and the loss 
rate j/2ir of the metamaterial is 2THz and 0.2 THz for 
(a) and (b), respectively. A height of 28 nm corresponds 
to a normalized distance of zo/Xo f=a 0.02, where we can 
reasonably expect the dipole approximation to hold (the 
spatial dependence will be shown below). 

Our material loss numbers (7/^0 = 0.01 — 0.001) are 
close to the state-of-the-art for metamaterials, but they 
are significantly greater than those used in some previ- 
ous waveguide studies, where enhanced PFs were demon- 
strated with 7/cjo w 10~ 10 - 10" 
all [23( 1 . For metamaterial applications, 
of merit (FOM), is FOM = -Re(n)/Im(n), with a larger 
FOM indicating a less lossy metamaterial. The current 
FOM for typical metamaterial is of order 100 at GHz 
frequencies and drops to 0.5 at optical frequencies (380 
THz) [25||. However, there are methods to improve these 
FOMs as they are not fundamental material properties; 
for example, Soukoulis et al. have suggested that the 
FOM can be improved by a factor of 5 at optical frequen- 
cies, and after optimizing their fishnet design, they have 
demonstrated that the FOM can be around 10 at 380 THz 
(cf. Fig. 5(c) in Ref. [24|). Recently, they also demon- 
strated a new design where the FOM is about 60 at 40 
THz [H] . For the proposed structure with 7/2tt = 2 THz, 
we have a FOM of 0.72 at the resonance frequency, and 
a maximum FOM of 26.25 at 220 THz, which is similar 
to the state-of-the art FOMs at optical frequencies. 

The PF (Eq. \T7\ is enhanced at the frequency of 
ui/2tt w 189 THz (w ujo/2it), and the enhancement for 
a z— polarized dipole is larger than that for an x or 
y— polarized dipole; part of the reason that the PF is 
larger for the z— polarized dipole is that it couples to the 
TM modes, which are more strongly influenced by ma- 
terial losses (through /i); we have verified that the TM 
and TE PFs approach one another as 7 — * 0, and this 
trend can partly be seen by comparing the TE and TM 
PFs in Fig. 4a and Fig. 4b. In the presence of the larger 
loss (7/2tt = 2 THz) the peak F z (TM PF) is about 240 
and the peak F x / y is about 120. When 7/27r is decreased 
to 0.2 THz, the PF increases significantly to T z = 720 
and T x / y = 350. The physical origin of these large PF 
enhancements comes from the slow energy velocities of 
the propagation modes. This can be seen from the dis- 
persion curves (Fig. 2 and Fig. 3), where upon close in- 
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FIG. 4: Purcell factor as a function of frequency for r y/2n = 
2 THz (a) and 7/2tt = 0.2 THz (b). The dot is located at 
20 = h/W = 28 nm (0.02 Ao) above the NIM slab. The blue 
curves correspond to a 2— polarized dipole, and the red curves 
correspond to an x/y— polarized dipole. The inset is the re- 
fractive index, n, vs frequency, where the dashed line cor- 
responds to the imaginary part and the solid line to the real 
part (which is negative throughout the entire frequency range 
shown) . 



spection, we realize that we are obtaining the odd mode 
resonances, as is expected from dipoles near the surface 
where the local field is larger near the NIM surface; for 
example, the resonance around 207 THz corresponds to 
the v g — > region of the TE3 mode (cf. Fig. 2(c) on the 
complex— uj band structure). In the complex—/? disper- 
sion curves (Fig. 3(b)), this same resonance is seen as 
the point where the two branches of the TE3 mode split 
apart near 207 THz ((3 ~ 7 fj,m). The series of peaks 
below 200 THz and above uj are due to the slow group 
velocity region of the various odd modes, which approach 
one another at too- 

In addition, we note that for the TE modes, there 
is a PF peak due to the plasmon resonance around 
223 THz that has also been highlighted elsewhere, e.g. 
Refs. [22|, [23|, H3] ■ What is particularly interesting, is that 
for reduced losses, this resonance splits into two (cf. Fig. 
4(b)); moreover, by inspection of the band structure (cf. 
Fig. 3(b)), the lower-lying peak of this pair is actually 
due to the bound propagation modes (e.g., TE5), which 
is only visible in the complex—/? band structure. For 
larger losses, these individual peaks cannot be resolved, 
and one must then assume that the broadened peak near 
223 THz for 7/27r = 2 THz, is due to a combination of 
the TE SPP and bound modes contained within the light 
lines. We emphasize that this resonance, which is below 
the SPP frequency, is not observable with the complex— to 
band structure, as discussed earlier; it is alo unique to the 
NIM structure. 

Next, the PF for different dot positions zq is investi- 
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FIG. 5: (a) The dependence of the maximum values in the PF 
on position zq for a loss factor of *y/2n = 2 THz.. The solid 
blue and red curves correspond to the PF at uq for TM and 
TE modes respectively and the dashed red curve corresponds 
to the TE SPP mode. The upper and lower dashed grey 
curves show the scaling of l/^o' 5 an d 1/zo, respectively. Even 
at the smallest distances from the NIM surface (~ 7nm), the 
electrostatic limit (where a scaling of l/«o would occur) has 
still not been reached at the chosen frequencies (see text), (b) 
The maximum PFs as a function of damping factor for the 
dot located at zq = h/20. The curve labelling is the same as 
in (a). All curves are obtained at a fixed frequency. 



gated. Similar to what happens near metal surfaces, as 
the dot is brought closer to the NIM slab surface, the 
PF will increase rapidly and formally diverge at the sur- 
face of the NIM slab. In fact, it is a straightforward 
exercise to show analytically that the electrostatic G at 
the surface of a half-space lossy structure has an imagi- 
nary part that diverges. One finds for the non-retarded 
terms (quasi-static approximation) [54|: G(r s ,r s ) = 
G°(r s , r s ) T G°(r s , -r s )(e nim - l)/(s n im + 1), where the 
minus and plus sign refer to TE and TM polarization 
respectively, G° is for free space, and e n im is the permit- 
tivity of the NIM medium. Because Re {G°(r s , r s )} di- 
verges, any amount of loss in s n i m , no matter how small, 
will lead to a divergent LDOS at the surface. Conse- 
quently, the quasi-static approximation will not work at 
the surface, and in general we should consider distances 
significantly larger than the emitter size if we are to em- 
ploy the dipole approximation. 

Figure OJb) shows the dependence of the values of the 
PF at three different resconance peaks on the position Zq. 
Because of the expected LDOS divergence at the surface, 
we only show the behavior down to distances of h/40 
(zq/Xq — 0.005); we expect the dipole approximation to 
work to distances of around h/10 (zo/Xq = 0.02). In 
obtaining these graphs we have fixed the frequency. The 
PF enhancements arc found to decrease as a function 
of distance, as expected, but large values can still be 
obtained at distances (~ OAh) or more. For example, 
the TM peak has a PF of 10 at a distance of 100 nm 
from the surface. For metal surfaces, the TM SPP mode 
PF scales as I/zq for small distances (e.g. see Ref. poj)). 



which we have also verified for our structure (see Fig. 
5(a)). This behavior is due to the electrostatic scaling 
of the evanescent contribution from the SPP, and for our 
structure, this scaling dominates for distances of around 
zq < 0.04 h. However, the scaling of the metamaterial 
modes is quite different: we obtain a scaling of around 
1/zp 1 - 5 for the PF peak at iv and also for the TE SPP 
mode peak. In the limit of fi = 1, we again recover the 
1/zq scaling for these modes. Also in the case of the 
metamaterial, significant PFs can still be achieved for 
much larger distances away from the surface, even for 
zq = h. The reason for this unexpected scaling is that 
the chosen resonance frequencies have not yet recovered 
the electrostatic limit, even for dipole distances as small 
as 7nm from the surface; if we choose frequencies away 
from the waveguide peaks, then we indeed get the the 
1/zq scaling from the NIM, as expected. 

The dependence of spontaneous emission on the damp- 
ing factor is plotted in Fig. [HJb), which shows that in- 
creasing the damping factor decreases the peak sponta- 
neous emission rate, while increasing the full width at 
half maximum (FWHM) of the PF resonance. How- 
ever, even in the presence of very large losses (e.g., 
7/27T =10 THz), we see that large PFs are still achiev- 
able. As expected, clearly there is also a large improve- 
ment in the PF enhancement if the nominal losses can be 
improved by an order of magnitude. 



B. Lamb shift and far-field spectrum of 
spontaneous emission 

The Lamb shift is another fundamental quantum effect 
whereby the vacuum interaction with a photon emitter 
can cause a frequency shift of the emitter [4l|. Cavity 
QED level shifts of atoms near a metallic surface have 
been shown to be significant as one approaches the sur- 
face [l9[ . For optimal coupling, usually these are studied 
at high frequencies (e.g., Tiui > 4cV), so as to couple to 
the TM SPP resonance; as a function of frequency, the 
level shift changes sign as we cross the resonance [l!| • For 
lower frequencxies, a 1/zq van der Waals scaling again 
occurs [551 ]. Given the complicated modal structures of 
NIM waveguides, it is not clear what the Lamb shifts will 
look like, and to the best of our knowledge, Lamb shifts 
have never been studied in the context of metamaterial 
waveguides. 

The QED frequency shift of the emitter can be di- 
rectly calculated via the real part of the scattered GF 
(Eq. (19)). Experimental dipole moments for semicon- 
ductor quantum dots vary from around 30 D (D=Debye) 
to 100 D [U, [HJ, so here we adopt a realistic dipole mo- 
ment of d = 50 D. The results for z— polarization and 
x/y— polarization are shown in Fig. HJa) and Fig. HJb), 
respectively. Our calculations indicate that the frequency 
for a dipole above the NIM slab will be significantly red- 
shifted relative to vacuum, with rich frequency oscilla- 
tions as one sweeps through the NIM waveguide reso- 
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FIG. 6: Frequency shift due to the inhomogeneous scatter- 
ing when r y/2n = 2 THz (solid curves) and y/2iv = 0.2 THz 
(dashed curves). The dipole strength is d — 50 D (see text), 
(a) The dot polarization is perpendicular to the slab surface 
(TM). (b) The dot polarization is parallel to the slab surface 
(TE). 



nances. The Lamb shift at the too resonance frequency 
for different loss factors 7 are identical, and are not 
zero; the nonzero shift at the main resonant frequency 
uiq is due to the asymmetry of the PF, namely the se- 
ries of waveguide resonances at the higher frequency 
end of loq- The frequency shift at luq for TM modes is 
6.3 GHz, and for TE modes is 3.5 GHz, and the differ- 
ence between them mainly comes from TM SPP modes 
at u) pe /\/2 w 490 THz. When the loss is reduced, the 
various modal contributions become more pronounced. 
It is worth highlighting that this frequency shift, which 
is completely unoptimized, is already comparable to some 
of the largest shifts reported for the real-index pho- 
tonic crystal environment, e.g., |&<j|/u; = 4 x 10~ 5 [53|. 
In normalized units, we obtain frequency shifts around 
I Su I /oj — 5 x 10~ 5 , over a wide frequency range. This ra- 
tio is even larger for smaller distances (and larger dipole 
moments), however one must watch that the dipole ap- 
proximation does not breakdown. 

We also remark that these NIM Lamb shift features are 
substantially different to those predicted in typical met- 
als. For example, we have calculated the Lamb shift from 
a half space of Aluminium and find that the Lamb shift 
in the same optical frequency regime is only -0.1 GHz, 
for an identical dipole and position. Although closer to 
the SPP resonance (which is at 2780 THz), much larger 
values (PF = 322 at zq = 7nm) can be achieved, the fre- 
quency dependence is relatively featureless, in contrast to 
that shown in Fig. 6 and the spontaneous emission is very 
strongly quenched in metals near the SSP resonance. 

Finally, we turn our attention to the spontaneous emis- 
sion radiation that can actually be measured. In the 
following, we use Eq. (21) to investigate the emitted 
spontaneous emission spectrum for two different exci- 
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FIG. 7: Detected spontaneous emission spectrum from a dot 
above the metamaterial (thick blue curve), compared with the 
emission from the same dot in free space (thin black curve). 
The dot location is zo = 28 nm above the surface and the 
detector position is 2Ao fa 2800 nm above the quantum dot. 
The dot polarization is perpendicular to the slab surface, (a) 
The excition frequency is is ujdi = 189.1 THz. (b) The excition 
frequency is u)d2 = 190.7 THz. 



ton frequencies, Ud = ujdi and LOd2 which are indicated 
in Fig. [HIa). The polarization of the quantum dot ex- 
citon is assumed to be along z, and the loss factor is 
7/27T = 2 THz. The frequency dependence of the emit- 
ted radiation is shown in Fig. [7] for the two different dot 
frequencies. For ujd = Udi, the shift at peak emissions 
is —6.3 GHz and the enhancement in spontaneous emis- 
sion rate is 240. For ojd = Wd2, the shift is —4.7 GHz and 
the enhancement in the spontaneous emission rate is 135. 
Usually, with a large dipole moment of d = 50 D and a 
spontaneous emission enhancement on the order of 100, 
the the photon-dot interaction will enter the strong cou- 
pling regime and the emitted spectrum will show a typical 
spectral doublet. However, since the emitted field away 
from the NIM is predominantly carried away by radiation 
modes, there the far field contains signatures of the QD 
coupling, even to near fields; however, because there is 
quenching in the far-field emission, strong coupling is not 
observable in the far field from the NIM waveguide-such 
effects could possibly be observed in the near field given 
suitable detection capabilities, e.g., from a scanning near 
field optical microscope. These features depend upon the 
properties of the GF propagator appearing the sponta- 
neous emission formula (Eq. (21)). Comparing with the 
free space emission spectrum, the integrated emission in 
Fig. 7(a) and Fig. 7(b) is 1.8 x 10~ 4 and 1.6 x 10" 3 , re- 
spectively. Thus, the predicted far field spectra, which 
obtain Purcell factor and Lamb shift signatures, should 
certainly be observable experimentally. 
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VI. CONCLUSION 

In summary, we have employed a rigorous medium- 
dependent theory and a stratified Green function tech- 
nique, to investigate the enhanced emission characteris- 
tics of a single photon emitter near the surface of the NIM 
slab waveguide in the optical frequency regime. The ori- 
gin of the predicted Purcell factor peaks is primarily due 
to slow light propagation modes which have been ana- 
lyzed by calculating the complex band structure of this 
waveguide. Correspondingly, we also predict a signifi- 
cant frequency (Lamb) shift of the single photon emitter, 
with rich features that stem from the waveguide mode 
characteristics. All of our predictions are based on a 



realistic metamaterial model that includes both mate- 
rial dispersion and loss and scales to any region of the 
electromagnetic spectrum; the role of material loss and 
dipole position has also been investigated in detail. It 
is further shown that the rich emission characteristics at 
the surface can act as a sensitive and non-perturbativc 
probe of below-light-line waveguide mode characteristics. 
These predicted medium-dependent QED effects are fun- 
damentally interesting and may find use for applications 
in quantum information science. 
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